function [dof_map, V, T, u1, u2, err] = ex2(n_round, FE_Type, FE_Order, Young, nu, method)
% function [dof_map, V,T, u1, u2, err] = ex2(n_round, FE_Type, FE_Order, Young, nu, method)
%
% The problem given in suri paper 
% 
%
%  Author: Dr. Xian-Liang Hu
%

Young = 250; 
nu = 0.25;
% lambda = 100;
% mu = 100;
R = 1/2;

if nargin == 0
    FE_Type = 'BB';
    FE_Order = 3;
end

Quad_Order = 13;


% generate initial mesh
[p, e, t] = initmesh('squareg', 'Hmax', 0.35, 'Hgrad', 1.99);
V = [p(1,:)' -p(2,:)']; T = t(1:3,:)';
[T, E, ET, TE] = build_fem_mesh(V, T);

% n_elem = size(T,1); n_dof_per_elem = (FE_Order + 1)*(FE_Order + 2)/2;

%%%%%%%%
%
% set the edge flags for boundary treatment
%
% edge_mark = zeros(size(E,1),1);  % default flag is inner flag, which is 0
bnd_flag = (ET(:,2)==0);
% edge_mark(bnd_flag) = -1;  % free boundary, homogeneous Neumann, which is -1

%% given the Dirichlet and Neumann boundary by hand, this is problem
%% dependent! 

%% this is only for test Dirichlet boundary!
bdr_Dirichlet = find(bnd_flag == 1);
% bdr_Nermann = [];


%%% this is for the linear FE case:
% n_dof = max(max(dof_map)); 
% u = zeros(n_dof,1); u(dof_map) = u1;
% trisurf(T,V(:,1),V(:,2),u);

%%% order d finite element solutions:
disp_d = 3; tri_temp = template_mesh_tri(disp_d);
[u_h, Tris1, Points1] = fe_solution_bb(V,T, u1, FE_Order, tri_temp, disp_d);
[v_h, Tris2, Points2] = fe_solution_bb(V,T, u2, FE_Order, tri_temp, disp_d);

% analytic solutions
[uu, vv] = fun_u(Points1(:,1),Points1(:,2), Young, nu, R);

% cal the L_inf error
erru = max(max(abs(u_h - uu)));
errv = max(max(abs(v_h - vv)));

% return the maximum
err = max(erru, errv);

% % visualization and print information
% subplot(1,2,1); 
trisurf(Tris1, Points1(:,1), Points1(:,2), u_h);  %plot error for u
% subplot(1,2,2); trisurf(Tris1, Points1(:,1), Points1(:,2), v_h - vv);  %plot error for v


fprintf('The inf norm or error(u) = %e,  error(v) = %e.\n',  erru, errv);

end


%%%%%%%%%%%%%%%%%%%%%%%%
%  analytic solution, also used as boundary function
% 
function [u_1, u_2] = fun_u(xx, yy, E, nu, R)
   r2 = (xx.*xx + yy.*yy);
   u_1 = r2/100 + (1 - 1/100)*R*R;
   u_2 = u_1;
end


%%%%%%%%%%%%%%%%%%%%%%%%
%  The right hand side function
% 
function [f_1, f_2] = fun_f(xx, yy, E, nu, R)
    [n_row, n_col] = size(xx);
    f_1 = -8*ones(n_row, n_col);
    f_2 = f_1;
end